# This files assess the robustness of Dahlgaard (2018) 
# by varying the bandwidth from three to 183 days

load("agg_day_year_all.rdata")
load("n_day_year_all.rdata")

## commented out for replication files 

#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("reshape2")
#  install.packages("dplyr")
#  install.packages("rdrobust")

library(xtable)
library(ggplot2)
library(dplyr)
library(rmeta)
library(rdrobust)
library(reshape2)

## Define treatment dummy and assign weights 
data_day_year <- 
  data_day_year %>%
  ungroup() %>%
  mutate(treated = days > 0)

#create large data set 
data <- left_join(data_day_year, data_n_year,
                  by = c("days", "year"))

data_close <- 
  data %>%
  mutate(voted     = round(par_vote*n),
         abstained = round((1-par_vote)*n),
         treated   = days > 0) %>% 
  select(year, days, treated, voted, abstained)

data_voted <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["voted"]])) %>%
  mutate(voted = 1)

data_abstained <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["abstained"]])) %>%
  mutate(voted = 0)

data_master <- 
  rbind(data_voted, data_abstained) 

## Loop over bandwidths

effect <- data_frame(yr2009 = rep(NA, 183),
                     yr2013 = rep(NA, 183),
                     yr2014 = rep(NA, 183),
                     yr2015 = rep(NA, 183),
                     pooled = rep(NA, 183),
                     place  = 1:183,
                     rd     = "Conventional RD")

se <- data_frame(yr2009 = rep(NA, 183),
                 yr2013 = rep(NA, 183),
                 yr2014 = rep(NA, 183),
                 yr2015 = rep(NA, 183),
                 pooled = rep(NA, 183),
                 place  = 1:183,
                 rd     = "Conventional RD")

effect_robust <- data_frame(yr2009 = rep(NA, 183),
                            yr2013 = rep(NA, 183),
                            yr2014 = rep(NA, 183),
                            yr2015 = rep(NA, 183),
                            pooled = rep(NA, 183),
                            place  = 1:183,
                            rd     = "Robust RD")

se_robust <- data_frame(yr2009 = rep(NA, 183),
                        yr2013 = rep(NA, 183),
                        yr2014 = rep(NA, 183),
                        yr2015 = rep(NA, 183),
                        pooled = rep(NA, 183),
                        place  = 1:183,
                        rd     = "Robust RD")

years <- c(2009, 2013, 2014, 2015)
for (i in 3:183){
  for(j in 1:4){
    sub_data <- 
      data_master %>% 
      filter(year == years[j])
    rd                  <- rdrobust(y = sub_data$voted, x = sub_data$days, h = i)
    effect[i, j]        <- rd$coef[1]
    se[i, j]            <- rd$se[1]
    effect_robust[i, j] <- rd$coef[3]
    se_robust[i, j]     <- rd$se[3]
  }
  effect[i, 5] <- meta.summaries(as.numeric(effect[i, 1:4]), as.numeric(se[i, 1:4]))$summary
  se[i, 5]     <- meta.summaries(as.numeric(effect[i, 1:4]), as.numeric(se[i, 1:4]))$se.summary
  effect_robust[i, 5] <- 
    meta.summaries(as.numeric(effect_robust[i, 1:4]), as.numeric(se_robust[i, 1:4]))$summary
  se_robust[i, 5]     <- 
    meta.summaries(as.numeric(effect_robust[i, 1:4]), as.numeric(se_robust[i, 1:4]))$se.summary
  
  if(i == 10*floor(i/10))
    print(i)
  
}

effect_joint <- 
  rbind(effect, effect_robust) %>% 
  melt(., id = c("place", "rd")) 

se_joint <- 
  rbind(se, se_robust) %>% 
  melt(., id = c("place", "rd")) %>%
  rename(se = value)

bw_data <- 
  data_frame(variable = c("yr2009",
                          "yr2013",
                          "yr2014",
                          "yr2015",
                          "pooled"),
             bw = round(c(with(data_master %>% filter(year == 2009),
                               rdbwselect(y = voted, x = days))$bws[1],
                          with(data_master %>% filter(year == 2013),
                               rdbwselect(y = voted, x = days))$bws[1],
                          with(data_master %>% filter(year == 2014),
                               rdbwselect(y = voted, x = days))$bws[1],
                          with(data_master %>% filter(year == 2015),
                               rdbwselect(y = voted, x = days))$bws[1],
                          NA)))

plot_data <- 
  left_join(effect_joint, se_joint,
            by = c("place", "rd", "variable")) %>%
  mutate(low = value + qnorm(0.025)*se,
         hi  = value + qnorm(0.975)*se) %>% 
  left_join(., bw_data, by = "variable")

labels <- c(yr2009 = "2009",
            yr2013 = "2013",
            yr2014 = "2014",
            yr2015 = "2015",
            pooled = "Pooled")
plot <- 
  ggplot(data = plot_data,
         aes(x = place, y = value)) + 
  geom_line() + 
  facet_grid(variable ~ rd,
             labeller = labeller(variable = labels)) +
  geom_line(aes(x = place, y = low), alpha = 0.5, linetype = 2) +
  geom_line(aes(x = place, y = hi), alpha = 0.5, linetype = 2) + 
  theme_classic() +
  xlab("Bandwidth in days") + 
  theme(axis.title = element_text(size = 14)) + 
  geom_vline(aes(xintercept = bw), alpha = 0.5, linetype = 2) + 
  geom_hline(yintercept = 0, alpha = 0.3, linetype = 3) + 
  scale_y_continuous("Effect in percentage points", limits = c(-0.15,0.15), breaks = seq(-0.1, 0.1, by = 0.05)) +
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14))

